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Two Reynolds-averaged Navier-Stokes codes were used to compute flow over the NASA 
Trapezoidal Wing at high lift conditions for the AIAA CFD High Lift Prediction 
Workshop, held in Chicago in June 2010. The unstructured-grid code FUN3D and the 
structured-grid code CFL3D were applied to several different grid systems. The effects 
of code, grid system, turbulence model, viscous term treatment, and brackets were stud- 
ied. The SST model on this configuration predicted lower lift than the Spalart-Allmaras 
model at high angles of attack; the Spalart-Allmaras model agreed better with experiment. 
Neglecting viscous cross-derivative terms caused poorer prediction in the wing tip vortex 
region. Output-based grid adaptation was applied to the unstructured-grid solutions. The 
adapted grids better resolved wake structures and reduced flap flow separation, which was 
also observed in uniform grid refinement studies. Limitations of the adaptation method as 
well as areas for future improvement were identified. 


I. Introduction 

The AIAA CFD High Lift Prediction Workshop (HiLiftPW-1) was held at the 28^^ Applied Aero- 
dynamics Summer Conference in Chicago, IL during June 2010.^ HiLiftPW-1 was patterned after the 
successful AIAA Drag Prediction Workshops.^ The objectives of HiLiftPW-1 were to assess the numerical 
prediction capability of current-generation CFD for a representative three dimensional (3D) configuration. 
The HiLiftPW-1 test cases were based on the Trapezoidal Wing configuration, which was used in a series of 
NASA Langley 14- by 22-Foot Wind Tunnel and NASA Ames 12-Foot Pressurized Wind Tunnel tests. 

The submissions for two NASA CFD codes are presented for simulations with several different grids as well 
as grid- adapted simulations based on these grids. The uniform grid refinement and adaptive grid refinement 
studies are primarily intended to ascertain the role of discretization error in the calculations. There are other 
sources of modeling error (e.g., geometry, wind tunnel corrections, boundary layer transition, and turbulence 
modeling) that limit the prediction of the wind tunnel measurements. Calculations with different turbulence 
models and brackets are also performed to better understand some of the remaining error sources. 

Chaffin and Pirzadeh^ detailed the effort required to manually specify grid resolution for accurate 3D 
high-lift computations. Mavriplis^ showed that discretization error still dominates other error sources (e.g., 
turbulence model distance function calculation, thin-layer viscous approximation, level of artihcial dissipa- 
tion) even with extremely hue grids that are computable with today’s supercomputers. To address these 
issues, adaptive grid methods are designed to automatically specify a resolution request that equidistributes 
and minimizes error estimates. 

The improved resolution request is commonly based on local error estimates. Uniformly reducing the 
errors associated with ah local-error sources of the flow may not be optimal from an engineering context, 
where calculating an output functional may be of greater concern. An alternative method is to estimate 
the error in the calculation of a specihed engineering output functional. Output error indicators utilize 
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the dual or adjoint solution of an output functional to account for the impact of local error as well as the 
transport of these local errors throughout the problem domain to improve the calculation of that output 
functional. 

The workshop summary and participant talks showed the greatest variation in the calculations for the 
flap and wing tip of the model. This paper is organized to begin with non-adaptive results to show trends 
in these sensitive flap and wing tip regions with uniform grid refinement, grid construction, geometry, and 
physical model. The adaptive grid results are then examined at these sensitive conditions to show common 
trends. 


II. CFL3D Flow Solver 

CFL3D^^ is a structured-grid upwind multi-zone CFD code that solves the generalized thin-layer or full 
Navier-Stokes equations. It can use point-matched, patched, or overset grids and employs local time-step 
scaling, grid sequencing and multigrid to accelerate convergence to steady state. A time-accurate mode is 
available, and the code can employ low-Mach number preconditioning for accuracy in computing low-speed 
steady-state flows. CFL3D is a cell-centered finite- volume method. It uses third-order upwind-biased spatial 
differencing on the convective and pressure terms, and second-order differencing on the viscous terms; it is 
globally second-order accurate. Roe’s flux difference-splitting method^^ is used to obtain fluxes at the cell 
faces. The solution is advanced in time with an implicit approximate factorization method. A wide variety of 
turbulence models are available in the code, including linear eddy-viscosity and nonlinear models. For each 
loosely coupled iteration, the mean flow equations are advanced in time with the eddy-viscosity fixed then 
the turbulence model is advanced in time with the mean flow solution fixed. Only the one-equation model 
of Spalart-Allmaras^^ (SA) and the two-equation shear stress transport (SST) model of Menter^^ were used 
for this study. 


III. FUN3D Flow Solver 

FUN3D^^”^^ is a finite- volume Navier-Stokes solver in which the flow variables are stored at the vertices 
or nodes of the mesh. FUN3D solves the equations on mixed element grids, including tetrahedra, pyramids, 
prisms, and hexahedra and also has a two-dimensional path for triangular/quadrilateral grids. At interfaces 
delimiting neighboring control volumes, the inviscid fluxes are computed using an approximate Riemann 
solver based on the values on either side of the interface. Several convective flux schemes are available in 
FUN3D. The most common scheme for subsonic and transonic flows is Roe’s flux difference splitting, which 
is used in the current study. For second-order accuracy, interface values are obtained by extrapolation of 
the control volume centroidal values, based on gradients computed at the mesh vertices using an unweighted 
least-squares technique. Several flux limiters are available in FUN3D to limit the reconstructed values when 
necessary. The most common used is that of Venkatakrishnan,^^ which is used in this study. For tetrahedral 
meshes, the full viscous fluxes are discretized using a finite- volume formulation in which the required velocity 
gradients on the dual faces are computed using the Green-Gauss theorem. On tetrahedral meshes this is 
equivalent to a Galerkin type approximation. For non-tetrahedral meshes, the same Green-Gauss approach 
can lead to odd-even decoupling. A pure edge-based approach can be used to circumvent the odd-even 
decoupling issue, but yields only approximate viscous terms. Thus for non-tetrahedral meshes, the edge- 
based gradients are combined with Green-Gauss gradients, which improves the h-ellipticity of the operator, 
and allows the complete viscous stresses to be evaluated. This formulation results in a discretization of 
the full Navier-Stokes equations. 

The solution at each time-step is updated with a backwards Euler time-differencing scheme. At each time 
step, the linear system of equations is approximately solved with either a multi-color point-implicit procedure 
or an implicit-line relaxation scheme. Local time-step scaling is employed to accelerate convergence to 
steady state. For turbulent flows, several models are available within FUN3D. In the current study, only the 
loosely coupled SA model^^ was used for FUN3D. 

IV. Output-Based Adaptation 

Output-based adaptation requires the solution of the discrete adjoint equations.^^’^^’^^ These linear 
adjoint equations are solved with a dual-consistent time-marching method. The iterative method is con- 
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stmcted such that an adjoint solution exactly dual to a direct differentiation approach is recovered at each 
time step, yielding a convergence rate which is asymptotically equivalent to that of the ffow equations. To 
improve adjoint iterative convergence, a generalized conjugate residual (GCR) scheme^^ is used with the 
standard multi-color point-implicit procedure as search directions. 

Venditti^^ described an output-based error estimation and adaptation scheme, which utilizes an embedded 
grid. This has been implemented in within a paralleff^ framework. The embedded grid required for 
the Venditti error estimator can be prohibitively large (eight times original grid) for 3D problems, which 
effectively limits the problem size of the output adaptive method. A single-grid error estimation procedure 
has been described and compared to the embedded-grid procedure for Euler problems. This is the first 
time the single- grid error estimate has been applied to turbulent simulations. 

Venditti^^ provided a procedure to calculate a new grid spacing request from the adaptive indicator and 
an user specified error tolerance. In this study, this error tolerance is set to half the total remaining error 
estimate at each adaptive iteration. The anisotropy of mesh elements is based on the Mach Hessian, where 
the element size in the smallest spacing direction is dictated by the output-based new grid spacing request. 
The anisotropic spacing request is specified by the metric M. 

Many large-scale simulations are resource limited; it may not be practical to compute a sufficient grid 
size required to reach a grid converged answer. A constant-complexity scaling^ can optionally be applied to 
M to allow the user to directly specify the size of the adapted grid. This size can be set to obtain an optimal 
grid that is computable with available resources. To perform this scaling, an estimate of the adapted grid 
size Nq is computed as 



where the integral is approximated by a sum over dual control volumes in the domain Q, each with a volume 
of Vn. The metric is scaled to produce a requested adapted grid size A^request by 

2 

( 2 ) 

This specified complexity adaptation is applied to the single-grid error estimation results. 

Parallel grid mechanics described by Park and Darmofal^^’^^ are used to modify the grid. The adaptation 
mechanics have only been implemented for the tetrahedral element type. The local elemental operators of 
node insertion, node movement, element swap, and element collapse are utilized to iteratively drive the edge 
lengths to unity in M. The grid adaptation mechanics suffer from robustness problems when performing 
highly anisotropic adaptation near curved boundaries. Also, the discretization in FUN3D is more accurate 
for semi-structured right-angle elements in the boundary layer. To address both of these issues, the near 
wall boundary layer grid is kept constant or frozen during the adaptation process. For this study, the grid 
is fixed within 0.14 inches of solid boundaries. This approach is not applicable to problems with initially 
under-resolved geometry or boundary layers, because that portion of the grid is held constant during the 
simulation. However, this approach is a powerful technique for problems with important off-body anisotropic 
features, i.e., vortices and wakes. Details of the frozen boundary layer adaptation technique are available in 
Andren and Park^^ and Park and Carlson. 

V. High-Lift Prediction Workshop Cases 

All of the HiLiftPW-1 workshop cases were based on a trapezoidal three-element high-lift wing (Trap 
wing) configuration with a simple single-slotted flap and slat. The model was developed in order to provide 
a database for CFD validation. This semi-span model has been tested both in the NASA Ames 12-Foot 
Pressurized Wind Tunnel (PWT)^^ and the NASA Langley 14- by 22-Foot Subsonic Wind Tunnel (SWT).^^ 
The SWT data was selected for the workshop. The photo of the model in the PWT, Fig. 1(a), provides a 
clear view of the body and brackets. The upper surface experimental pressure taps are shown in Fig. 1(b), 
with the locations examined in this study highlighted. Although multiple slat and flap settings were tested, 
the workshop test cases focused primarily on a flap landing configuration (Configuration 1). The nominal slat 
angle for Configuration 1 was 30 degrees, and the nominal flap angle was 25 degrees. The reference geometric 
quantities for the Trap wing are summarized in Table 1. Details regarding the experimental geometry, data 
and repeatability can be found in Slotnick et al.^^ or on the HiLiftPW-1 website. While the wind tunnel 
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model was tested as a semi-span model, all the workshop test cases were specified to be run in “free-air” 
using a computational symmetry plane at the center line of the geometry. Similarly, no wind tunnel walls 
or model support systems were to be included except as specified below in Case 3. The experiment did not 
have specified transition on the model. However in the workshop, the boundary layer was to be modeled as 
“fully turbulent “ in the CFD simulations i.e. turbulence is modeled in all areas of computational domain. 
All computations were performed at a freestream Mach number of Mo© =0.2 and freestream chord Reynolds 
number of Rcc = 4.3 x 10^. 

The first test case (Case 1) for the workshop was a grid convergence study for Configuration 1 at two 
angles of attack: a = 13° and 28°. A family of three grids was required for the study with consistent levels 
of refinement between each grid: coarse, medium, and fine. An additional extra-fine mesh solution was 
optional. The second test case (Case 2) for the workshop was a flap deflection study. The flap deflection 
study required an angle-of-attack sweep for the Trap wing with two different flap deflection angles: 25° 
(Configuration 1) and 30° (Configuration 8). For the flap deflection study, computations were performed 
on a set of medium-sized grids. Angle of attack computations were performed at 6°, 13°, 21°, 28°, 32°, 34° 
and 37°. The optional slat and flap support effects study (Case 3) was a comparison of Case 1 clean wing 
results on the medium grid size with results from a grid including the slat and flap brackets. Rumsey et al.^^ 
provides a summary of all of the workshop submissions. 



(a) Photo of model in PWT (b) Location of upper surface pressure taps 


Figure 1. Trap wing wind tunnel model. 


Table 1. Reference geometry for the Trap Wing. 


Mean Aerodynamic Chord (c^e/) 

39.634 inch 

Wing Reference Area 

22.028 ft2 

Wing Span 

85.054 inch 

Aspect Ratio (AR) 

4.561 

Quarter Chord Sweep angle 

29.97 deg 
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VI. CFL3D and FUN3D Non-Adaptive Results 

Several complete sets of runs were made using CFL3D and FUN3D on committee-supplied grids, and 
were submitted to the HiLiftPW-1 workshop. In addition, a few new runs (not a complete submission) on 
a grid with brackets created by participants from JAXA were made using FUN3D for this paper. All of 
the non-adaptive runs are summarized in Table 2. The FUN3D results on tetrahedral grids NCTetA were 
made as a preliminary step toward computing adaptive results, as described in the next section. Details 
about the specific grids are not provided here, but documentation can be found in Slotnick et al.^^ or on 
the HiLiftPW-1 website. 


Table 2. Summary of non-adaptive runs. 


Code 

Full Grid Name 

Abbreviated 
Grid Name 

Turbulence 

Model 

Brackets? 

Viscous Term 
Treatment 

CFL3D 

Str-OnetoOne-A 

OnetoOneA 

SST 

No 

thin-layer* 

CFL3D 

Str-OnetoOne-A 

OnetoOneA 

SA 

No 

thin-layer* 

CFL3D 

Str-OnetoOne-B 

OnetoOneB 

SST 

No 

thin-layer*-l- 

CFL3D 

Str-OnetoOne-A 

OnetoOneA 

SA 

No 

full N-S 

FUN3D 

U nst- Tet- N o decent er ed- A 

NCTetA 

SA 

No 

full N-S 

FUN3D 

Unst-Mixed-FromTet- 

Nodecentered-A 

Mix(NCTetA) 

SA 

No 

full N-S 

FUN3D 

Unst-Mixed-Nodecentered-C 

NCMixC 

SA 

Yes 

full N-S 


* = full N-S with viscous cross-derivatives neglected 
-h = results not included in this paper 

Many workshop participants noted a sensitivity of their solutions to the initial conditions of the com- 
putation, especially at the higher angles of attack. The current authors also found that initializing the 
computation with free-stream conditions could result in a solution with stall characteristics at lower angles 
of attack than a solution that was initialized from a previously converged computation at a lower angle of 
attack. Multiple solutions sometimes occurred at the higher angles near stall. The results presented in this 
paper represent the solution sets that produced higher stall angles. 

At the workshop, the CFL3D results presented using SST on the OnetoOneB medium grid stalled at a 
lower angle of attack than the solutions using the OnetoOneA medium grid. It was recently noted that all 
of the OnetoOneB computations were initialized with free-stream conditions even at the higher angles of 
attack near stall. Since the effect of initial conditions on the OnetoOneB grid solutions were not quantified, 
OnetoOneB results will not be included in this paper. 

VI. A. Grid Convergence Study (Case 1) Forces and Moments 

The effects of grid refinement were explored for Configuration 1 at o = 13° and 28°. Fig. 2 shows Cl, Cl>, 
and Cm as a function of where N represents the number of grid points (for a node-centered scheme) 

or grid cells (for a cell-centered scheme). On sufficiently refined grids that lie in the asymptotic range of grid 
convergence, a spatially second-order method should produce a straight line when plotted in this manner. 
As seen in the figure, the behaviors were not always linear. The SA results with different codes and grids 
tended to come together to within about 0.05 in Cl, 0.005 in Cl>, and 0.015 in Cm as the grid was refined 
for a = 13°. The spread was much larger for a = 28°. Experimental values of the forces and moments 
interpolated to these flow conditions are shown as a filled dot on the ^-axis of the individual figures. For the 
most part, grid refinement improved agreement with experiment. 

For the FUN3D solutions, the tetrahedral grid, NCTetA, was more sensitive to grid density than the 
mixed-element grid base on the same nodes, Mix(NCTetA). For both Mix(NCTetA) and NCTetA, the 
coarse grids yielded results that appear to be out of the asymptotic range. Use of a different family of 
mixed-element grids (NCMixC) tended to yield slightly lower Cl values than Mix(NCTetA) despite the fact 
that the NCMixC grids had a globally higher number of nodes. 
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(a) a = 13°, Cl 


(b) Q = 28°, Cl 
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--0-- 

CFL3D, SA, OnetoOneA 
CFL3D, SST, OnetoOneA 
CFL3D, SA, OnetoOneA FNS 
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■ FUN3D,SA, Mix(NCTetA) 
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experiment 

X 

• 
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2E-05 3E-05 
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(c) Q = 13°, Cd 


(d) a = 28°, Cd 




(e) a = 13°, Cm 


(f) a = 28°, Cm 


Figure 2. Grid convergence, Configuration 1. 
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For the CFL3D solutions, the use of full Navier-Stokes produced higher lift, higher drag, and more 
negative moment than the approximation that neglected viscous cross-derivatives. Thus the inclusion of 
the viscous cross-derivatives consistently improved the correlation with the experimental force and moment 
results. Other effects of neglecting viscous cross- derivative terms will be shown in the next section. The 
grid convergence study also showed that the CFL3D results using SST converged toward values that were 
significantly lower in lift, lower in drag, and less negative in moment than the SA results. Further comparisons 
of the SA and SST turbulence model results will also be shown in the next section. 

VI. B. Modeling Sensitivities 

Figure 3 shows the lift curves for all of the CFL3D computations for Configuration 1 on the medium grids. 
Similar to what was noticed at the workshop by other HiLiftPW-1 participants, computations using the SST 
model generally were found to yield lower lift levels (in poorer agreement with experiment) than computations 
using the SA model particularly near Cl, max- Typically, SST produced more flap separation than SA. An 
example is shown in Fig. 4, which shows surface-restricted streamlines for Configuration 1 at a = 13° using 
CFL3D. Only SA results will be shown in the remainder of the paper, because the SA model results correlated 
much better with experiment than the SST results. 



Figure 3. Lift curves for CFL3D computations, Configuration 1, medium grids. 



(a) SA (b) SST 

Figure 4. Surface-restricted streamlines, CFL3D, a = 13°, Configuration 1, medium OnetoOneA grid. 

Additionally, Fig. 3 shows some differences between computations that used full Navier-Stokes (FNS) and 
those that neglected viscous cross-derivative terms: the full Navier-Stokes solutions yielded higher lift levels 
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near Cl, max- Here, we explore a key reason for the differences. Fig. 5 shows surface pressure coefficients at 
two locations for Configuration 1 at a = 28°. In the figure, “TLNS” indicates the run with viscous cross- 
derivative terms neglected. At 98% span location, the TLNS solution yielded dramatically different surface 
pressures in much poorer agreement with experiment than full Navier- Stokes. The pressure coefficients at 
the flap forward span location (shown earlier as “flapfwdspan” in Fig. 1(b)) show that the difference occurred 
primarily in the tip region. Although FNS did not capture the proper Cp levels near the wing tip, TLNS 
was much worse. As will be shown below for the FNS case, there was only a relatively small influence of 
grid refinement near the wing tip, but the trend was in the correct direction. 




(a) 98% span station on main element 


(b) Flap forward span station 


Figure 5. 


Surface pressure coefficients, Configuration 


1, a = 28°, CFL3D, SA, fine grid. 


Fig. 6 shows vorticity contours at 4 different x-planes as a way to visualize the tip vortex development 
(note that small gaps seen in the geometry are due to a plotting package anomaly and are not real). In the 
TLNS case where the viscous cross-derivative terms were neglected, the tip vortex was more diffuse than the 
case where full Navier-Stokes was employed. With a more compact tip vortex, the resulting FNS pressures 
yielded greater suction in the tip region, in better agreement with experimental data. This difference 
demonstrates that the use of full Navier-Stokes is very important for this configuration; thin- layer- type 
approximations should be avoided. 



(a) Viscous cross-derivatives neglected (b) FNS 


Figure 6. Vorticity contours showing tip vortex development, CFL3D, Configuration 1, a = 28°, SA, medium OnetoOneA 
grid. 
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VI. C. Slat and Flap Support Effects Study (Case 3) 

The effect of brackets was investigated by using the NCMixC set of grids from JAXA. Only two angles of 
attack {a = 13° and 28°) were run on medium grids. Computations with brackets yielded lower lift by 
ACl = 0.0138 at o = 13° and by ACl = 0.0519 at o = 28° which moves the computational results away 
from the experimental values. Other workshop participants showed a similar effect for Case 3. Results at 
the flap forward span station and at the 50% flap span station are shown in Figs. 7(a) and (b), respectively. 
With brackets modeled, the pressure coefficients along the flap showed a “scalloping” effect in Fig. 7(a), 
similar to experimental results. At the 50% span station on the flap, shown in Fig. 7(b), including brackets 
yielded a slight improvement in the agreement with experimental surface pressure coefficients. (Note other 
workshop participants saw a greater improvement in the agreement with experimental surface pressures at 
this span station.) 




(a) Flap forward span station (b) 50% span location on flap 

Figure 7. Surface pressure coefficients showing the effect of brackets, Configuration 1, a = 28°, FUN3D, NCMixC grids, 
SA. 


The influence of the brackets on the a = 28° flowfleld behind the wing is shown in Fig. 8. There was 
some local roll-up of vorticity due to the presence of the brackets. However even for the clean wing, the 
wake structures tend to diffuse quickly when compared to the structured grid results in Fig. 6, which have 
a C-mesh topology with inherent wake resolution. The effect of mesh adaptation on the wake structure will 
be shown in the next section. 




(a) No brackets (b) Brackets included 

Figure 8. Vorticity contours at x = 60” (6/ac/c), 70” (g^reen), 80” (6Zne), and 90” (red) constant planes, showing the effect of 
brackets. Configuration 1, a = 28°, FUN3D, NCMixC grids, SA. 
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The effects of both grids and brackets can be visualized by plotting surface restricted streamlines from 
FUN3D, shown in Fig. 9 for a = 13°. At this angle of attack, the flap trailing edge had separation along 
much of the span (recall earlier CFL3D results in Fig. 4). Here, the effect of different grids can be seen by 
comparing Figs. 9(a) and (b). The use of the NCMixC grid produced somewhat less flap separation than 
use of the Mix(NCTetA) grid. The reason for this difference is not known. Comparing Figs. 9(b) and (c), it 
can be seen that the presence of the brackets caused significant modulation of the flap separation line. The 
effect of flap angle is studied in the next section. 




(a) Mix(NCTetA) grid, No brackets 


(b) NCMixC grid, No brackets 



(c) NCMixC grid. Brackets included 


Figure 9. Surface restricted streamlines showing the effect of grid and brackets, Configuration 1, a = 13°, FUN3D, SA. 


VI. D. Flap Deflection Study (Case 2) 

A subset of the runs will be focused on for the remaining portion of this paper based on experience gained 
from studying the modeling sensitivities. Figure 10 shows lift curves, drag polars, and moment polars for 
CFL3D on the OnetoOneA grid, FUN3D on the NCTetA grid, and FUN3D on the Mix(NCTetA) grid. 
All methods shown used the full Navier-Stokes formulation and SA turbulence model. For Case 2, angle 
of attack sweeps for Configuration 1 and Configuration 8 were performed on the medium grids. The left 
column shows Configuration 1 and the right column shows Configuration 8. All three sets of computational 
results were very similar, with the largest differences in the moment coefficient. The lift and drag coefficients 
were generally predicted slightly low compared to experiment, and the moment coefficient was generally 
predicted low in magnitude compared to experiment. Experiment measurements are only available at and 
below a = 28° for Configuration 8 Recalling Fig. 2, refining the grid tended to bring computed results closer 
to experimental levels. For Configuration 1, the computations consistently predicted stall between a = 34° 
and a = 37° with the slat stalling first followed by the wing. For Configuration 8, the CFL3D results showed 
a similar stall behavior, but the FUN3D results on the NCTetA grid stalled at a lower angle. The FUN3D 
results on the Mix(NCTetA) grid did not indicate the same type of stall mechanism for Configuration 8 as 
for Configuration 1. For the FUN3D computation on Mix(NCTetA), the wing stalled first. 
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a 


(a) Configuration 1, lift curve 


(b) Configuration 8, lift curve 




Co 


(c) Configuration 1, drag polar 


(d) Configuration 8, drag polar 




(e) Configuration 1, moment polar 


(f) Configuration 8, moment polar 


Figure 10. Force and moment curves, medium grids. 
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Figure 11 shows the delta force and moment of Configuration 1 minus Configuration 8. The delta lift, 
drag, and pitching moment curves indicate that the effect of flap deflection was slightly under-predicted 
by the computations at the lower angles of attack and over-predicted at the higher angles. The CFL3D 
results were closer to the experimental trends than FUN3D at the higher angles of attack. For the FUN3D 
results, the inconsistent stall angle prediction between Configuration 1 and Configuration 8 led to an erratic 
prediction of the flap deflection effect near the stall angle. 




a a 


(a) Delta lift curve 


(b) Delta drag curve 



a 


(c) Delta moment curve 

Figure 11. Delta force and moment curves, medium grids. 


VI. E. Grid Convergence Study (Case 1) Surface Pressures 

Comparing the surface pressures from the grid convergence study illustrate the areas of the flow where the 
code-to-code and grid-to-grid are significant. The grid convergence of FUN3D and CFL3D surface pressures 
are shown for the 85% span station of the flap (Configuration 1 and a = 13°) Fig. 12 FUN3D exhibited 
greater grid sensitivity than CFL3D. The CFL3D variation was localized at the peak suction of the flap. 
The FUN3D refinement showed increased suction on the forward portion of flap and a change in the flap 
trailing edge separation characteristic. Although not shown, use of FUN3D with all tetrahedral NCTetA 
grids yielded even greater grid sensitivity. However, on the finest grid the results overall were very similar. 
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(a) FUN3D on Mix(NCTetA) grids (b) CFL3D on OnetoOneA grids, FNS 

Figure 12. Grid convergence of surface pressure coefficients, Configuration 1, a = 13°, SA, 85% span station on flap. 


The flap forward span station (Configuration 1 and a = 28°) is is shown in Fig. 13. As seen in Fig. 1(b), 
the flap forward span station is near quarter chord of the flap upper surface. Both FUN3D and CFL3D have 
an improvement in the prediction of the tip vortex with grid refinement. FUN3D exhibited greater variation 
due to grid refinement in the inner portion of the wing. 




Y 


(a) FUN3D on Mix(NCTetA) grids 


(b) CFL3D on OnetoOneA grids, FNS 


Figure 13. Grid convergence of surface pressure coefficients. Configuration 1, a = 28°, SA, flap forward span station. 


Fig. 14 compares four different combinations of grid and code on fine grids for Configuration 1, a = 13°, 
at the 85% span station on all elements. Results were almost indistinguishable on the slat and main elements, 
but FUN3D produced slightly less upper surface suction and a greater amount of trailing edge separation on 
the flap than CFL3D. The effect of adaptive mesh refinement on this latter characteristic will be explored 
in the next section. 

Similarly, Fig. 15 compares the Configuration 1, a = 28°, fine grid results for the 85% span station on all 
elements. At the 85% span station, all four sets of results on fine grids showed only minor differences with 
the greatest variation in the flap upper surface peak pressure. In Fig. 15(c), all of the CFD predictions had 
a lower pressure on the flap lower surface than the measurements. The reason for this discrepancy is not 
known. 


13 of 24 


American Institute of Aeronautics and Astronautics 



Q. 

O 



50 


55 


60 

X 


65 


70 



X 

(a) Slat (b) Main (c) Flap 

Figure 14. Surface pressure coefficients, Configuration 1, a = 13°, SA, 85% span station on all elements, fine grids. 
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Figure 15. Surface pressure coefficients. Configuration 1, a = 28°, SA, 85% span station on all elements, fine grids. 
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A similar comparison of fine results at the flap forward span station is also shown in Fig. 16. There is 
a slight variation in pressure coefficient results across the flap span with the tip and root regions showing 
the greatest variation. All of the computations predict a lower peak-suction value near the tip region than 
measured in the experiment. 



Y 


Figure 16. Surface pressure coefficients, Configuration 1, a = 28°, SA, flap forward span station, fine grids. 


VII. FUN3D Adaptive Results 

Each of the adapted cases use a converged solution on the original grid at that angle of attack. These 
original grid solutions were computed with the aid of angle of attack sequencing; each angle of attack solution 
was restarted from a converged lower angle of attack solution on the same original grid. Attempts to start 
the high angle of attack cases from freestream initial conditions (not shown) were problematic. These high-a 
cases started from freestream either failed to converged or appeared to converge to a condition with massive 
separation and lower coefficient of lift than the a-sequenced cases. Adaptation on these massively separated 
solutions did not reattach the flow or reach a significantly higher coefficient of lift. Angle of attack sequencing 
was a practice employed by many workshop participants. 

Coefficient of lift is plotted as a function of coefficient of drag for the workshop-provided coarse and 
medium node-based tetrahedral Configuration 1 grids in Fig. 17(a). The original grid is show as open 
symbols at the workshop-specified angles of attack. Each of the original grids was adapted and the results 
are shown as filled symbols, Fig. 17(a). The adapted coarse grids used the embedded-grid error estimate 
and the adapted medium grids used the single-grid error estimate. The single-grid error estimate adaptation 
included constant complexity metric scaling, Eq. (2). Adaptation with the embedded- and single-grid error 
estimates is compared for the medium grid at a = 13° in the following sections. The adapted coarse grid 
results approximated the original medium grid results at a = 28° and below. The maximum lift was under- 
predicted by the original and adapted coarse grid computations. When looking at the entire polar, there 
was little difference between the original and adapted medium grid results. Drag and lift were lower for the 
high-(a pre-stall adapted grids, but were on the same curve as the original grids. 

VII. A. Coarse Grid Embedded- Grid Error Estimate Adaptation 

A zoom of the lift-drag polar at o = 13° is shown in Fig. 17(b). The coarse (3.6M nodes), medium (IIM 
nodes), and fine (32M nodes) original grids are shown as the open circles. The coarse and medium original 
grids were used for adaptation with the embedded-grid error estimate (filled square) and the single-grid 
error estimate (open squares). The series of adapted grids are connected by solid lines, forming a trajectory. 
The script driving the adaptive processes begins by restarting the flow solver, so the first point of each 
adaptive trajectory is at a slightly different location than the uniform refinement points due to different 
iterative convergence levels. The embedded- grid error estimate adaption based on the coarse original grid 
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Figure 17. Lift as a function of drag, Configuration 1. 


is approaching the original medium grid point. The final adapted coarse grid has IIM nodes. The medium 
grid adaptations, which are shown in more detail in the next section, also trend toward the fine grid point. 
Output-based adaptation is formulated to address discretization error and not modeling error, so it should 
tend toward the uniformly-refined grid answer and not the experimental measurements when modeling errors 
are present. 

The original and final adapted volume tetrahedral grid for the coarse grid adaptation at a = 13° is 
sliced at the 85% span station. Fig. 18. The intersection of a plane and the tetrahedral elements is shown 
for the original. Fig. 18(a), and adapted. Fig. 18(b), grids. The resolution of the slat, main element, and 
flap stagnation stream lines has been improved. The upper surface of the main element and wakes have 
been resolved with anisotropic elements. The i^-velocity at the 85% span station is shown for the original. 
Fig. 18(c), and adapted. Fig. 18(d), grids. The velocity deficit in the wake is better defined for the adapted 
grid. The adapted slat and flap gap regions are detailed in the next figures. 

The slat gap region dX a = 13° is sliced at the 85% span station. Fig. 19. The adapted grid. Fig. 19(b), 
is more refined in the slat stagnation, gap, and wake regions than the original grid. Fig. 19(a). The region of 
frozen boundary layer grid can be seen close to the wing in Fig. 19(b). The slat gap jet has a more defined 
boundary in the adapted grid. Fig. 19(d). The slat wake velocity deficit is also better defined and propagated 
a longer distance over the main element in the adapted grid. Fig. 19(d), than the original grid. Fig. 19(c). 
The flow has greater acceleration above the slat wake as compared to the original coarse grid solution. 

The flap gap region dX a = 13° is sliced at the 85% span station. Fig. 20. The adapted grid. Fig. 20(b), 
is more refined near the flap stagnation region, gap, upper surface, and wake. The adapted grid is also more 
refined near the main element lower surface, cove, and upper surface. The region of frozen boundary layer 
grid can also be seen close to the wing in Fig. 20(b). The peak velocity in the flap gap jet increases for 
the adaptive grid as well as the region above the main element wake. The region of strongly reversed flow 
above the aft portion of the flap (large negative values of rt- velocity) is decreased in the adaptive grid. These 
trends are also indicated by the grid convergence of the surface pressure coefficients for the Mix(NCTetA) 
grids. Fig. 12(a). With grid refinement, the peak negative pressure coefficients increases indicating a higher 
velocity at x=ll and the region of separated flow is decreased at x=76. 

VII. B. Medium Grid Single- and Embedded-Grid Error Estimate Adaptations 

Details of the medium grid adaptive trajectories are provided in Fig. 21(a) as a closer zoom of Fig. 17(b). 
The adaptive trajectories for the single- and embedded-grid error estimate were similar. Due to the resources 
required for the embedded-grid error estimate, this is the only case performed to compare the single- and 
embedded- grid error estimates. The final embedded- grid adaptation is 17M nodes and the final single grid 
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(c) Original coarse grid n- velocity. (d) Final adapted coarse grid n- velocity. 


Figure 18. Slice of volume grid at the 85% span station for the original and adapted coarse grid, Configuration 1, 
a = 13°. 
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(c) Original coarse grid n- velocity. (d) Final adapted coarse grid n- velocity. 


Figure 19. Slice of volume grid at the 85% span station for the slat of the original and adapted coarse grid, Configuration 
1, a = 13°. 
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(c) Original coarse grid n- velocity. (d) Final adapted coarse grid n- velocity. 


Figure 20. Slice of volume grid at the 85% span station for the flap of the original adapted coarse grid, Conflguration 
1, a = 13°. 
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adaptation is 23M nodes. 

The simulation adapted with the single grid estimate became unsteady as seen in the iterative convergence 
plot, Fig. 21(b). The spikes in residual were due to high frequency errors in the interpolated solution of the 
adapted grid that quickly damped out. The range of coefficient of lift and coefficient of drag are the same 
for Fig. 21(a) and Fig. 21(b). The lift and drag are oscillating through a significant range of this zoomed-in 
figure. This case was also run in a time-accurate mode (not shown) but this did not reduce the range of 
the oscillations. This type of small scale unsteadiness can lead to incorrect sensitivities^^ because a variable, 
pseudo-time formulation was used for the flow and adjoint solver, which assumes a steady solution. This 
problem is exacerbated by the adaptive scheme increasing grid resolution and reducing dissipation induced 
by insufficient grid resolution. For all cases, the iterative convergence of the SA turbulence model stalled on 
adapted grids after two adaptations. The source of this stalled iterative convergence is still being investigated. 
The convergence histories of lift and drag at other angles of attack (not shown) reached steady-state. 
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(b) Convergence history of adaptation with single-grid error 
estimate. 


Figure 21. Medium grid adaptation, Configuration 1, a = 13°. 


The coefficient of pressure on the flap of the original medium grid and medium grid adapted with the 
single-grid error estimate are shown in Fig. 22 for the 85% span station at a = 13°. The adapted grid data is 
a snap shot of an unsteady simulation. Fig. 21(b). The adapted grid better approximates the experimental 
coefficient of pressure measurements than the original grid. The indication of separation on the initial medium 
grid at X = 74.5 is reduced with adaptation. The trend of reduced separation with grid refinement are also 
indicated by the grid convergence of the surface pressure coefficients for the Mix(NCTetA) grids. Fig. 12(a). 
There is less flap separation indicated by the adapted grid than the fine tetrahedral and mixed-element grids. 
Fig. 14(c). 

The surface pressure coefficient is shown in Fig. 23 for the original and adapted grids at a = 28°. This 
is the same condition as in Fig. 5 and Fig. 13. As seen in Fig. 23(a), the upper surface pressure decreased 
at the 98% span station of the main element with grid refinement, showing improved comparison to the 
experimental measurements and FNS structured grid results. Fig. 5(a). The negative coefficient of pressure 
peak moved outboard with a slightly more negative value for the adapted grid. Fig. 23(b). The location 
of the peak was well predicted, but overall levels were less negative than the measurements. This is the 
opposite of the trend indicated by the unstructured and structured uniform refinement studies. Fig. 13. The 
difference in trends with grid refinement may be due to the frozen geometry and near wall boundary layer 
grid of the adaptive results. The small peak near y = —6 due to the side of body vortex is better resolved. 
It is also clearly seen in the structured grid results. Fig. 13(b). 

Vorticity magnitude contours at x-constant volume slices are shown in Fig. 24 for the original and adapted 
grids at a = 28°. The structured grids were sliced at the same locations shown in Fig. 6. These same slices 
were shown from a different perspective in Fig. 8 for unstructured grids with and without brackets. The 
vorticity rapidly dissipated downstream for the original medium grid. Fig. 24(a). The vorticity in the 
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Figure 22. Medium grid flap coefficient of pressure at 85% span station, Configuration 1, a = 13°. 




(a) Main element at the 98% span station 


(b) Flap at forward spanwise location 


Figure 23. Surface pressure coefficients, Configuration 1, a = 28°, original and adapted grids. 
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boundary layer, wake, and vortices was better defined in the adapted grid. Fig. 24(b). The vorticity in the 
slat wake was better defined for the adapted grid, and the tip and side-of-body vortices were more clearly 
defined. The original grid slat wake was almost indistinguishable from the main element and flap boundary 
layers, whereas the adapted grid created improved resolution of the wakes, similar to the C-mesh topology 
of the structured grid. Fig. 6. 



(a) Original medium grid. (b) Adapted medium grid. 


Figure 24. Vorticity contours showing tip vortex development, FUN3D, Configuration 1, a = 28°, SA, medium grid. 


This study has highlighted the need for improvements in the current implementation of the output-based 
adaptation process. The frozen boundary layer limitation of the adaptation mechanics may be preventing 
the calculation of the fine grid result starting with a coarse initial grid. Fig. 17(b). The a = 13° case became 
unsteady as the grid was refined. Fig. 21(b), which indicates the need for grid adaptation schemes compatible 
with time-accurate simulation. Monitoring the iterative convergence and adaptation process for the large 
number of cases required to compute the angle of attack sweep can be improved. The turbulence model 
failed to converge on the adapted grids, which will require further examination. 

VIII. Conclusions 

The contributions of two Reynolds-averaged Navier-Stokes codes - one structured and one unstructured - 
to the 1®^ AIAA CFD High Lift Prediction Workshop were described. Two types of submissions were made: 
(a) a set of standard runs on coarse, medium, and fine grids, exploring the effects of grid, code, turbulence 
model, and viscous term treatment; and (b) a set of output-based grid-adapted solutions, exploring the 
influence of automated grid refinement. These submissions follow the current state-of-the-art for CFD 
applied to high-lift applications. The latter submissions represent active development at the leading edge of 
CFD capabilities. 

Among the standard runs, there were several different combinations of grid, code, turbulence model, 
and viscous term treatment. The majority of these were official workshop submissions, while one set of 
computations on the grids provided by JAXA (including brackets) was only performed for this paper. The 
SST turbulence model was found to produce more separation and lower lift than the SA model. The latter 
model’s results agreed better with experiment. Two independently-generated structured grids were found 
to yield different behaviors near Cl, max when using the SST model, although the particular aspects of the 
grids that might have caused the difference were not discovered. 

The effect of grid refinement was assessed for a = 13° and 28°. Among all the runs that used the SA 
model, the forces and moments did not approach a converged infinite-grid solution. In other words, the use 
of different grids and different codes yielded a non-negligible spread in results even as the grids were refined. 
The spread was greater at the higher angle of attack. A tetrahedral grid exhibited somewhat greater grid 
sensitivity than a mixed element version of the same grid (for which boundary layer tetrahedra were merged 
into prisms). The use of the approximation neglecting viscous cross derivative terms yielded lower lift than 
full Navier-Stokes. A key difference was in the prediction of the wing tip vortex; the thin-layer type of 
approximation yielded a more diffuse vortex and hence not as much suction force over the wing tip upper 
surface region. The use of full Navier-Stokes yielded results that agreed much better with experiment. The 
unstructured code tended to exhibit more trailing edge flap separation than the structured code at a = 13°, 
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but the particular unstructured grid used also had an effect on this ffow feature. The effect of brackets was 
explored in a qualitative sense. Including brackets yielded lower lift levels moving the results away from the 
measurements, but produced some features that bore greater similarity to experiment, such as a spanwise 
“scalloping” of the pressures over the ffap. 

The output-adaptive process was driven with both embedded-grid and single-grid error estimates. The 
single-grid error estimation procedure also included the option of adapting at a constant number of control 
volumes, which yields an optimized grid for a given computational resource. The near- wall regions of the 
grids were frozen due to limitations in the grid adaptation mechanics. The embedded-grid error estimate was 
applied to the entire set of workshop angles of attack with the workshop-provided coarse grid. The single- 
grid procedure was applied to the same conditions with the medium workshop grid. Both error estimation 
procedures were compared at a = 13° for the medium grid. Grid adaptation showed improved resolution of 
velocity, surface pressures, and vorticity magnitude. Changes in these quantities indicated a better resolution 
of tip vortices, side-of-body vortices, gap jets, and wakes for the adapted unstructured grids. The extent of 
boundary layer separation on the upper surface of the ffap is inferred by ffow velocity and surface pressures. 
The trend of reduced separation with adaptive grid refinement compares favorably with the uniform grid 
refinement studies. Limitations discovered in the current output-adaptive approach were detailed to guide 
further development of the method. 
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